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We show that Monte Carlo sampling of the Feynman diagrammatic series (DiagMC) can be used for tackling 
hard fermionic quantum many-body problems in the thermodynamic limit by presenting accurate results for the 
repulsive Hubbard model in the correlated Fermi liquid regime. Sampling Feynman's diagrammatic series for 
the single-particle self-energy we can study moderate values of the on-site repulsion (U /t ~ 4) and temperatures 
down to T/t = 1/40. We compare our results with high temperature series expansion and with single-site and 
cluster dynamical mean-field theory. 
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Advancing first-principle simulations of fermionic many- 
particle systems is notoriously difficult due to an exponential 
computational complexity jVL One of the most famous exam- 
ples is the Hubbard model |2|,|3|] whose phase diagram remains 
highly controversial, 
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where 6L creates a fermion with spin projection a =f, J, on 
site i, fiia = c' ia c ja , (...) restricts summation to neighboring 
lattice sites, and t, U and /u are the hopping amplitude, the 
on-site repulsion, and the chemical potential respectively. 

The Hubbard model is of great technological and scien- 
tific importance since it is believed by many to be the cen- 
tral model for high-temperature superconductors [3]. More- 
over, it can nowadays be realized with ultracold atoms in op- 
tical lattices where Mott physics JH101 has been observed 
recently. The prospect of using cold atomic gases as build- 
ing blocks of quantum simulators is thus becoming a real- 
ity. However, claiming the equation of state, either numeri- 
cally or experimentally, can only be substantiated on the basis 
of unbiased tools. Not only numerical methods need to be 
tested against known unbiased results, but also quantum em- 
ulators need to be validated by benchmarking them at ever 
lower temperatures. In the bosonic case this can be done with 
high precision ]7[], but for fermions it is still an open ques- 
tion what method can be used to accurately describe the low- 
temperature regime of the model. Precise numerical solutions 
of the Hubbard model can only be obtained in the particle- 
hole symmetric case = 1/2, where the sign problem is 
absent in quantum Monte Carlo (QMC) simulations. In all 
other cases circumventing exponential scaling necessarily in- 
volves approximations with systematic errors that are hard to 
assess and control. 

In general terms, Diagrammatic Monte Carlo (DiagMC) 18D 



is a numerical approach based on stochastic sampling of Feyn- 
man diagrams JsJGlil] to increasingly higher orders in the cou- 
pling constant — rather than evaluating integrals over internal 
variables for each diagram, one samples their momenta, imag- 
inary times and expansion orders stochastically. So far, Di- 
agMC algorithms have been successfully used for obtaining 
unbiased solutions for polaron problems (for the latest work, 
see Ref. ifTTll ). 



Nevertheless, the crucial question remained unclear of 
whether the DiagMC approach could be applied to true many- 
body systems in the thermodynamic limit at physically inter- 
esting temperatures much smaller than the Fermi energy. The 
convergence of a diagrammatic series for the many-body sys- 
tem is not guaranteed — in many cases the series has zero con- 
vergence radius, which can be revealed by Dyson's system 
collapse argument IU2I1 . Fortunately, the argument does not 
apply to the Hubbard model. However, even then the con- 
vergence is not guaranteed a priori due to possible unphysi- 
cal singularities in the parameter space, as, e.g., is the case 
with HTSE lll3U20ll . which for the Hubbard model diverges at 
T ~f. 

An algorithm for the many-body problem, based on a Di- 
agMC sampling of the self-energy, was developed recently by 
some of us 11411 . In this Letter, we find that this algorithm 
can produce controlled and accurate results (with reliable er- 
ror bars) in the correlated Fermi liquid regime of the Hub- 
bard model at moderate interactions in any spatial dimension. 
Comparison with single-site [15] and cluster [16-3 dynam- 
ical mean-field theory in the dynamical cluster approximation 
(DCA) variant in 2D and with high-temperature series expan- 



sions (HTSE) IU3L |2O0 reveals an impressive agreement, thus 
proving theoretical control over the Fermi liquid regime of the 
model at moderate values of U /t. Our simulations demon- 
strate that the diagrammatic expansion for the Hubbard model 
does converge (at least for moderate interactions) down to 
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temperatures much smaller than the Fermi energy i28ll . 

Despite the inherent sign problem, our scheme is capable of 
reaching diagram orders sufficient to claim the convergence 
in the thermodynamic limit, the feasibility of which for a 
many-body fermionic system was not clear a priori. The abil- 
ity to calculate higher-order diagrams is also important, e.g., 
in describing decay processes in QCD (see, e.g. 12111 ). and 
might prove relevant for solving Dyson-Schwinger equations 
in QED and QCD 12211 . where there is a growing interest in 
numerical perturbative expansions I23I1 . 

The crucial difference between the thermodynamic-limit 
DiagMC and other Monte Carlo approaches — dealing with 
a finite-size cluster of volume V in the (d + 1) -dimensional 
space — is that in DiagMC the complexity of the sign prob- 
lem coming from sign-alternating diagrammatic expressions 
is linked to the expansion order of the Feynman perturbation 
series and not to the large volume V required for the ther- 
modynamic limit extrapolation. Away from criticality and 
for moderate interaction strength, only a few diagram orders 
are needed until convergence of the diagrammatic series is 
reached (see Fig. Q]). This greatly reduces the computational 
cost, the major part of which is nothing but ensuring that 
higher order diagrams cancel each other and thus do not con- 
tribute to the answer. 

To demonstrate the power of DiagMC we will compare 
it to other state-of-the-art first-principle approaches: high- 
temperature series expansion and dynamical mean field theory 
(DMFT) in the single-site and dynamic cluster approximation 
(DCA) implementation HUSl. DMFT and its cluster ex- 
tensions have so far arguably been the method of choice for 
studies of the Hubbard model, producing a pha se diagram in 
close agreement with that of the cuprates 11911 . Single-site 
DMFT treats the correlated many-body system as a quantum 
impurity problem where each site can be viewed as an impu- 
rity embedded self-consistently in a reservoir produced by the 
other sites. It goes beyond conventional mean-field theory by 
taking all temporal correlations into account, but neglects the 
momentum dependence of the self-energy E. Short-range cor- 
relations can be taken into account by extending the size of 
the impurity, leading to various cluster schemes. 

HTSE is a method where the partition function Z = 
Tr exp(— $H) is expanded in (3f around the atomic limit (t = 0) 
and all the terms are summed up exactly up to a maximum ex- 
pansion order. This expansion is analytic and yields therefore 
precise thermodynamics of the Hubbard model, provided the 
series converges at the maximum expansion order considered. 

We consider the two cases of moderately (20%) and heavily 
(40%) doped regimes and moderate values of the interaction 
strength (U/t = 4) in two and three dimensions. The conver- 
gence of DiagMC is checked by performing simulations with 
ever increasing expansion order cut-off until the extrapolation 
to infinite order can be done with confidence. Although the 
sign problem is manageable for moderate interactions, simu- 
lations do get exponentially harder when the expansion order 
is increased. Fortunately, we can extract the converged re- 
sult before the sign problem becomes too severe. In case of 




FIG. 1: (color online) Energy density per spin component for the 
2D Hubbard model as a function of inverse expansion order cut-off 
N*. This example is showing data for U/t = 4 and T/t = 0.025. The 
density is shown in the inset. We observe convergence around order 
4, and the extrapolation to the infinite order can be done reliably. 
Error bars are smaller than the symbol size. 



the Hubbard model we typically see convergence of the se- 
ries around fourth order for a moderate interaction U/t — 4 
and sufficiently far away from half filling (see Fig. [T). All the 
results shown in this work are obtained by observing a clear 
plateau, from which we deduce the infinite-order answer with 
confidence. 

Let us first focus on the two-dimensional heavily hole- 
doped system with (n a ) w 0.3. Figure|2]shows the energy den- 
sity per spin component E a /tV as a function of temperature at 
a constant chemical potential p/t = —0.15 and on-site repul- 
sion U/t = 4. The density per spin component « a is shown 
in the inset. The lowest temperature in the plot, (3f = 40, is 
about two orders of magnitude below the Fermi energy e F . 
We clearly observe the Fermi liquid behavior, and the Fermi 
liquid parameters can be extracted by fitting the data to the 
expected the T 2 dependence (solid line in the figure) 



E(T)-E(0) = (p F + p' F e F ) 



K 2 T 2 



n{T)-n{0) = p' F 



, K 2 T 2 



(2) 

where p F and p^- are the Fermi-surface density of states 
and its derivative respectively. We find p F = 0.06(4), p' F = 
0.095(20). Note that in this highly overdoped regime the 
single-site DMFT results coincide with DiagMC data over the 
entire temperature range. 

Next, we consider the same system on the particle-doped 
side at a smaller doping of 20% doping ((« c w 0.6)). Figure[3] 
shows particle density at p/t = 3.1. In contrast to the previous 
case, we clearly see a systematic deviation of the single-site 
DMFT data from the DiagMC results for temperatures T <t. 
Upon increasing the cluster size the DCA data progressively 
approach the DiagMC points. Still, even 4x4 clusters are not 
sufficient to attain the accuracy of DiagMC in this parameter 
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FIG. 2: (color online) Energy density per spin component and den- 
sity (see inset) as a function of temperature for the 2D Hubbard 
model at fixed interaction strength U/t = 4 and chemical potential 
\x/t = —0. 15. Within error bars, DMFT (blue circles) coincides with 
the DiagMC (red squares) results. The Fermi liquid fit, see Eq.l(2j, is 
shown by solid lines. 
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FIG. 3: (color online) Comparison between the DiagMC and DCA 
methods: Density as a function of temperature for the 2D Hubbard 
model at U /t = 4 and p/t = 3.1. For temperatures T < t a systematic 
deviation of DCA data from DiagMC points is seen. The deviation is 
diminishing with increasing the size of the DCA cluster, the four-site 
results show deviations from single site and larger clusters because 
of the establishment of local short-range order 12511 . 



regime and an extrapolation in the cluster size is needed to 
reproduce the DiagMC results. 

The difference between single-site DMFT and DiagMC re- 
sults is reflected in the momentum dependence of the self- 
energy. By construction, E is momentum independent in 
single-site DMFT. In Fig. |4] we consider E(p,£) at the lowest 
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FIG. 4: (color online) Momentum-dependence of the self-energy 
at the Matsubara frequency % = 7t/f5 along the cut (0,0) — (71,0) — 
(71, jc) — (0, 0) in the first Brillouin zone for the Hubbard model with 
the parameters U/t = 4, fi/t = 3.1 and T/t = 0.4. DiagMC (squares) 
includes the full momentum-dependence, whereas single-site DMFT 
(solid lines) is momentum independent. The results of DCA calcula- 
tions are plotted by open symbols for clusters of size 4, 8, 16, and 32. 
Note that a good agreement between DiagMC and DCA is reached 
with 32-site clusters. The mean-field contribution (the Hartree term 
Una ~ 2.3f) was subtracted to magnify fine details. The arrows indi- 
cate the position of the Fermi momentum pf. 



Matsubara frequency t, = 7c/p along the cut (0,0) — (Jt,0) — 
(jt,Jt) — (0,0), in the (p x ,Py) plane. Cluster DMFT simula- 
tions do include momentum dependence of the self-energy ap- 
proximately. As expected, the DCA results get systematically 
closer to the DiagMC curves when increasing the cluster size. 
For a cluster of 32 sites, DCA shows a momentum dependence 
which is in good agreement with DiagMC. 

We now turn to the comparison with high-temperature se- 
ries expansion up to tenth order in (3f. To demonstrate that 
the DiagMC is insensitive to the spatial dimension, we will 
perform this comparison in three dimensions (3D). Figure [5] 
shows the energy and density dependence on temperature for 
the 3D Hubbard model at U/t = 4. In this particular exam- 
ple we keep the renormalized chemical potential pi = /j — Un a 
fixed. Within DiagMC, the variable pi is more convenient than 
the physical chemical potential p because the Hartree term 
Un c automatically accounts for tadpole diagrams, thus ex- 
cluding them explicitly from the simulation. Qualitatively, the 
data is similar to the 2D case. Without Pade approximation the 
bare HTSE series give an unbiased answer for temperatures 
only above T ~ 2t, while the degenerate Fermi liquid behav- 
ior develops only at T < t . Only then the energy and density 
display the characteristic quadratic dependence on tempera- 
ture. The fitted values of the density of single-particle states 
at the Fermi energy and its derivative are given in the caption 
of Fig. 

The data reported here can serve as established benchmarks 
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FIG. 5 : (color online) Energy density per spin component E a /tV and 
density n a (inset) as a function of temperature for the 3D Hubbard 
model at U /t = 4 and p' = p — Un a = 1.5. DiagMC results are shown 
by red squares. The HTSE data are shown by lines: second order 
in pf (green dashed line), eighth order (blue dotted-dashed line), and 
tenth order (red dotted line). At low temperature, the Fermi liquid be- 
havior is fitted to the DiagMC results (black solid line). The extracted 
Fermi liquid parameters are p(£f) = 0.26(4), p'(e/?) = -0.025(2). 



for present experiments on ultracold atoms in optical lattices 
and state-of-the-art numeric techniques. This work constitutes 
a proof of principle that DiagMC is a reliable method for deal- 
ing with hard many-body problems. 

To map out the phase diagram of the system one will need 
to sample diagrams for static response functions and observe 
the leading instabilities. In the degenerate Fermi liquid regime 
observed here, the development of these instabilities is con- 
trolled by the Fermi liquid theory, which allows one to reliably 
predict T c without actually simulating exponentially low tem- 
peratures. Further algorithmic progress can be made by using 
Dyson and/or Bethe-Salpeter equations to arrive at the self- 
consistent skeleton-diagrams description 112411 . which would 
allow one to consider larger interaction strengths and, pos- 
sibly, reach the critical points of the model. Going below 
the critical point might be possible by introducing anomalous 
propagators. 
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